Thermodynamic relations in a driven lattice gas: numerical experiments 
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We explore thermodynamic relations in non-equilibrium steady states with numerical experiments 
on a driven lattice gas. After operationally defining the pressure and chemical potential in the driven 
lattice gas, we confirm numerically the validity of the integrability condition (the Maxwell relation) 
for the two quantities whose values differ from those for an equilibrium system. This implies that a 
free energy function can be constructed for the non-equilibrium steady state that we consider. We 
also investigate a fluctuation relation associated with this free energy function. Our result suggests 
that the compressibility can be expressed in terms of density fluctuations even in non-equilibrium 
steady states. 

PACS numbers: 05.70.Ln,05.40.-a 



A rich variety of non-equilibrium phenomena have 
been successfully described by phenomenological evolu- 
tion equations. However, the microscopic foundation 
of such equations has not yet been established, except 
for systems near equilibrium. Even for non-equilibrium 
steady states (NESS) realized in simple systems, such 
as those involving only heat conduction and shear flow, 
appropriate statistical measures of microscopic configu- 
rations are not known. Recalling that equilibrium statis- 
tical mechanics was constructed with the aid of thermo- 
dynamics, we expect that checking the validity of ther- 
modynamic relations in NESS is an important step in 
constructing a theory of non-equilibrium statistical me- 
chanics. 

Non-equilibrium lattice gases are simple mathemati- 
cal models which have been useful in the elucidation of 
universal properties of NESS 0. Topics studied with 
such models include nonequilibrium phase transitions 2] , 
long-range spatial correlations 0, fluctuation theorems 
0, and non-local large deviation functionals p|, as well 
as mathematical foundations of nonequilibrium statis- 
tical mechanics It is thus expected that the non- 
equilibrium lattice gases provide good models for the ex- 
ploration of thermodynamic relations. 

There have been some proposals of an extended ther- 
modynamic framework applicable to NESS [3,|8|]- In one 
such study, Sasa and Tasaki start from operational def- 
initions of the pressure p and chemical potential fi, and 
they derive from these a quantitative relation which can 
be tested experimentally [g- Because the Maxwell rela- 
tion for p and fj, plays an essential role in the derivation 
of this relation, we are led to study the same Maxwell 
relation in the case of a driven lattice gas. 

In this Letter, we present numerical results that con- 
firm the validity of the Maxwell relation for p and fj,, 
which we define operationally for the system we study. 
As we explain below, the Maxwell relation provides an 
integrability condition for p and /z, and this yields a free 
energy function extended to the NESS that we consider. 
The existence of this free energy function leads us to 
believe that there is an associated fluctuation relation. 
Indeed, our numerical experiments suggest that the com- 



pressibility can be expressed in terms of density fluctua- 
tions even in certain non-equilibrium systems. 

Model: Let <Ji be an occupation variable defined on 
each site i = (i x , i v ) in a two-dimensional square lattice 
{(i x ,iy)\0 < ix < L+ 1,0 < i y < M + 1}. The variable 
o~i is 1 when the i-th site is occupied, and is if it is 
empty. We assume periodic boundary conditions in the 
x-direction ( i.e. Oi — o~j when when i = j + [L, 0)), and 
no flux boundary conditions in the y-dircction (i.e. <jj = 
when i y = 0, M+l). The array of all occupation variables 
{crj} is denoted as a and called the "configuration" . 

We study a driven lattice gas with the Hamiltonian 



H(a) = -) j o- l <7j - Ey^j x <jj, 

(hj) > 



(1) 



where represents a nearest neighbor pair and E is an 
external force 0] • The time evolution of a is described by 
the following rule: At each time step, choose randomly 
a nearest neighbor pair and exchange the values of 

<7j and o~j with the probability c(i,j; a) given by 



c(i,j;o-) 



l+aqp(p[H(o*t) -H{a)})' 



(2) 



where a v is the configuration obtained from a through 
this exchange and (3 is the inverse temperature. This ex- 
change probability is called the heat bath method and is 
one of the most standard update rules satisfying the local 
detailed balance condition, which condition is regarded 
to be natural in physical systems. The particle number 

= J2i a i is conserved throughout the time evolution. 
In this study, we fix f3 = 0.5, in which case the system is 
far from critical in the high temperature region. 

Pressure: In the equilibrium case (E = 0), the pres- 
sure is calculated by using equilibrium statistical mechan- 
ics. However, because we do not know the proper statis- 
tical measure in NESS, in this case we should define the 
pressure in an operational manner. Although the pres- 
sure is usually defined as the normal force exerted on a 
unit area of a surface, there is no quantity corresponding 
to the force in a lattice gas. We thus define the pressure 
in terms of the quasi-static work required to change the 
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FIG. 1: Pressure as a function of the density averaged at 
the center of the system. The circular and triangular symbols 
correspond to the cases E — and E — 10, respectively. Here, 
M = L = 32, Aw = 0.4 x 10" 5 , w m = 30 and to = 16000 
(msc). The averaged values for 100 samples are plotted. The 
statistical error bars are smaller than the symbols. The inset 
displays the L dependence of the pressure with N/L 2 = 0.5 
and L/M = 1. The pressure seems to converge to a definite 
value in the thermodynamic limit. 



system size. In order to allow for the calculation of this 
quantity in our lattice gas |2j , we add a wall potential to 
the Hamiltonian H(a) in the form 



,(<t) = H(a) + ^2 < J i w - 



(3) 



-M 



Because of the boundary conditions we impose, no parti- 
cles cannot exist at sites with i y = M + 1, and therefore 
the system size in the y direction is M. Now, according 
to © , as w is increased, it becomes increasingly unlikely 
for particles to exist at sites with i y = M . When w be- 
comes sufficiently large, the average occupation number 
for sites with i y = M can be considered zero. We de- 
note the value of w beyond which this is the case as w m . 
Hence, in the process that w changes from to w m , the 
effective system size in the y direction changes from M 
to M — 1. The quasi-static work performed to the sys- 
tem through this process is interpreted as the pressure 
multiplied by L. That is, the pressure p is written as 



p = lim — 

W m — >OG Jj 



dwJ^P 



NM 



dH w {a) 
dw 



(4) 



where Pjy M is the steady state distribution for a given 
value of w. It is easily proved that this definition of the 
pressure is equivalent to statistical mechanical formula in 
equilibrium. 

In the numerical experiments, values of p are obtained 
in the following way. Starting from random initial con- 
ditions, we carry out the time evolution with w = 
for a sufficiently long time, say t$. In this way, we 



obtain a steady state with w = 0. Next, we increase 
the value of w by a quantity Aw per Monte Carlo step 
per site (mcs) until w reaches w m . Then, noting that 
dH w {a)/dw = _ -m tTj , we measure the value of 
Sri =m ai as n (0 a t ti me *j where zero of the time is 
defined as the point at which w starts to increase. In one 
process from w — to w = w m , we calculate 



P 



Aw 



(5) 



4=1 



where t m is the time at which w — w m . Then, determin- 
ing carefully how the statistical distribution of p depends 
on w m and Aw, we estimate the value of p in the limit 
w m -> oo, Aw ^ 10]. 

In Fig. ^ we display an example of measured values 
of the pressure for densities in systems with E = and 
E = 10. It is important to note that our analysis is 
not restricted to systems near equilibrium. Indeed, the 
system with E = 10 is close to the strong field limit, 
in which the particle current is saturated to a constant 
value, and the equation of state for E — 10 clearly de- 
viates from the equilibrium one. This difference shows 
that the statistical distribution in the y direction differs 
from the equilibrium one. Also, the pressure becomes an 
intensive variable in the thermodynamic limit, as seen in 
the inset of Fig. 

Chemical potential: The chemical potential is mea- 
sured by placing a particle reservoir in contact with the 
system in the direction transversal to the driving field. 
We first assume that the chemical potential of the reser- 
voir, Hk(T, pr), is known. We also assume that there 
exists a chemical potential /i in this NESS that is a func- 
tion of T, p and E. From the definition of the chem- 
ical potential, we should have n(T,p,E) — nn(T, pr). 
Then, using p,(T, 1 — p, E) = pn(T, 1 — pr), which holds 
due to the particle-hole symmetry, the equality /i(T, p — 
1/2, E) = /iR,(T, pr, = 1/2) is derived. Thus, by measur- 
ing dp J (T,p,E)/dp for all values of p without contacting 
a particle reservoir, we can determine the form of the 
function /x(T, p, E). 

In order to measure dp,(T, p, E)/dp numerically, we 
add the one-body potential term J^. to the Hamil- 
tonian H(a), where cf>i — A(j> for i y > M/2 + 1, and 
4>i = for i y < M/2. We then measure the density 
profile along the y direction, in which there are two flat 
regions, 1 < i y < M/2 and M/2 + 1 < i y < M. We 
denote the density in the region 1 <C i y -C M/2 as p\ 
and the density in the region M/2 + 1 C i 9 C M as p2- 
When A(f) is sufficiently small, the chemical potentials of 
the two regions are equal by the definition of the chem- 
ical potential. By taking into account the shift in the 
potential energy, this condition can be written as 



t i(T,p 1 ,E)=p J (T,p 2 ,E)+Acf>. 

We thus obtain 

dp(T,p,E) 



dp 



= hm — — , 

A0^O Ap 



(6) 
(7) 
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FIG. 2: dp(T,p,E)/dp versus pdp(T, p, E)/dp for several 
values of p and -E. The triangle, square, star, plus, and circle 
corresponds to (p, E)=(0.5, 10), (0.4,10), (0.3,10), (0.5,3.0) 
and (0.5,0.0), respectively. Here, L = 32, and M = 32 for 
dp(T,p,E)/dp and M = 64 for pdp(T, p, E)/dp. The statis- 
tical error bars are smaller than the symbols. The dotted line 
corresponds to the Maxwell relation. 



FIG. 3: dtL£ as a function of l/M. Here, L = M = 32. The 
statistical error bars are smaller than the symbols. The dotted 
line is obtained as the best fit line of the form ai/M + b in 
the region l/M £ [0.4, 0.8], and d* is evaluated as the value of 
b. The inset displays the L dependence of d* with the values 
N/L 2 = 0.5 and L/M — 1 fixed, d* seems to converge to a 
definite value in the thermodynamic limit. 



where Ap = pi — P2 and p — (pi + p2)/2. Measuring 
Ap and p for several values of A0, we can evaluate the 
right-hand side of J7J) jlfl • 

Free energy: We begin by conjecturing that the rela- 
tion 



dp(T,p,E) dp(T,p,E) 
= P- 



dp 



dp 



(8) 



which holds in equilibrium, holds also in our NESS. This 
is equivalent to the Maxwell relation, because p and p are 
numerically confirmed to be intensive. Here, we estimate 
the value of dp/ dp (at p = 0.5, for example) by calculat- 
ing values of p for several values of p in a small interval 
around p = 0.5 10]. The results summarized in Fig. [21 
suggest the validity of the equality (JHJ). If indeed this 
relation does hold, its implication is significant, because 
when wc define the quantity F(T, M, N, E) as 



F(T, M, N, E) = N, [T, E) —MLp {t^e}, 

(9) 

p and p become 

V 
/' 



N 



1 dF(T, M, N, E) 
dM : 
dF(T, M, N, E) 



(10) 
(11) 



That is, F(T, M, N, E) can be regarded as the free en- 
ergy extended to the present NESS. Using the single 
function F(T, M, JV, E), we can derive various thermo- 
dynamic relations, including the Clapeyron law, under 
non-equilibrium conditions. 



Fluctuation relation: Consider the strip Qi — 
{(ix,*v)\l <ix< L,M/2-£/2- 1 < i y < M/2 + 1/2}, 
and define the density variable pi = X^en, °"-</|^|> wn i cn 
is coarse-grained over the strip fij. Let £ be the corre- 
lation length of density fluctuations in the y-direction. 
We define the free energy density f(p) for fixed (T,E) 
by F(T,M,N,E) = Mf(N/M) in the thermodynamic 
limit. It is then conjectured that the probability distri- 
bution of the density pi with ^<f<M can be written 
as 



P(p e = p) ~ exp(-/^[/(p) - /(p)]), 



(12) 



where p is the thermodynamic density. In the equilibrium 
case, such a form can be derived from a fundamental 
principle of statistical mechanics. Although there is no 
general proof of this form in the case of NESS, l|12|) seems 
plausible in the present case, because we have been able 
to define a free energy [l^. 

Recently, for some non-equilibrium lattice gases, the 
large deviation functionals of density fluctuations have 
been rigorously derived in non-local forms p|. These 
nonlocal forms are related to long-range correlations that 
exist generically in NESS [ll|. Here we shall not discuss 
this important issue further, and simply state our numer- 
ical finding that the scaling form 112fl has been clearly 
observed provided that one examines the density fluctu- 
ation in the strip f2f. 

If (|12fl is valid (at least locally), the fluctuation relation 







dp _ 1 
~d~ P ~ U({p-pf) 



(13) 



can be derived. This relation is known to be valid for 
describing fluctuations about equilibrium states, but it 
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FIG. 4: d* versus /3 x dp/dp for several values of (p,E). 
The triangle, square, star, plus, and circle correspond to 
(p,E) = (0.5,10), (0.4,10), (0.3,10), (0.5,3.0) and (0.5,0.0), 
respectively. Here, L = 32, and M = 32 for d„ and M = 64 
for P~ 1 dp/dp. The statistical error bars are smaller than 
the symbols. The dotted line corresponds to the fluctuation 
relation. 



is not known whether it is valid in NESS. To study this 
point, we investigate 1(13(1 numerically. 
We first note the asymptotic form 



di = (Pi) - (plY 



1 M-t 



(14) 



Li M 

for £ < I < M (see Fig. [3J, where is defined. Ac- 



cording to this, Lt((p — p) 2 ) in (fH!|> should correspond 
to <i*. The values of and f3~ 1 dp/dn are plotted in Fig 
. This result suggests the validity of 1(13(1 • In addition, 
combining 113(1 with (jSJ, the relation we obtain between 
the compressibility and density fluctuations is the same 
as that existing in equilibrium. 

We remark here that the asymptotic form 1(14(1 can 
be understood by considering the following simple sit- 
uation. Consider n random variables Xi (1 < i < n), 
with the conservation constraint X)j=i x % = 0- The sta- 
tistical properties of xi are given by E{xj) = and 
E(xiXj) — Sij, where E(x) represents the expected value 
of the random variable x. Let Xk be the partial sum of 
k elements randomly chosen from the set {x{\. Then, 
the probability of Xk with large k and large n — k is 
given by the Gaussian distribution with E(Xk) = and 
E(Xl) = k{n-k)/n 13]. This supports the form of (H). 

Discussion: We have presented new relations ob- 
tained using numerical experiments on a driven lattice 
gas. We hope that the present study can be extended to 
a wider variety of systems. Also, finding a connection be- 
tween the entropy, which is defined as -~(dF/dT)M,N,E 
in our formulation, and the Shannon entropy would be 
important in future construction of a statistical mechan- 
ical theory of NESS [l|. 

The authors acknowledge H. Tasaki for stimulating dis- 
cussions on steady state thermodynamics and driven lat- 
tice gases. They also thank Y. Oono and A. Shimizu for 
fruitful discussions on NESS. This work was supported 
by grants from the Ministry of Education, Science, Sports 
and Culture of Japan, No. 14654064. 



[*] Electronic address: hayashi@jiro.cu-tokyo.ac.jp, 
sasa@jiro.cu-tokyo.ac.jp 

[1] B. Schmittman and R. K. P Zia Phase transitions and 
Critical Phenomena Volume 17, Statistical Mechanics of 
Driven Diffusive Systems (Academic Press, 1995) ; J. 
Marro and R. Dickman Nonequilibrium Phase Tran- 
sitions in Lattice Models (Cambridge University Press, 
1999). 

[2] S. Katz, J. L. Lebowitz, H. Spohn, J. Stat. Phys. 34, 497, 
(1984). 

[3] P. Garrido, J. L. Lebowitz, C. Maes and H. Spohn, Phys. 

Rev. A 42, 1954 (1990). 
[4] J. L. Lebowitz and H. Spohn, J. Stat. Phys. 95 333, 

(1999). 

[5] B. Derrida, J. L. Lebowitz and E. R. Speer, Phys. Rev. 

Lett. 89, 030601 (2002) and references therein. 
[6] H. Spohn, Large Scale Dynamics of Interacting Particles 

, (Springer- Verlag, 1991). 
[7] Y. Oono and M. Paniconi, Prog. Theor. Phys. Suppl. 

130, 29 (1998), and references there in. 
[8] S. Sasa and H. Tasaki, cond-mat/0108365 
[9] We consider the pressure only in the direction transver- 



sal to the driving field. In general, the pressure tensor 
becomes anisotropic in NESS. However, we do not know 
how to measure the pressure in the direction parallel to 
the driving field. 

[10] K. Hayashi and S. Sasa, in preparation. 

[11] J.R. Dorfman, T.R. Kirkpatrick, and J.V. Sengers, Ann. 
Rev. Phys. Chem. 45, 213 (1994) and references therein. 

[12] The large deviation functional of density fluctuations in 
the model we consider turned out to be shape-dependent 
[F. J. Alexander and G. L. Eyink, Phys. Rev. E 57, 6229 
(1998).] . We have clarified that the large deviation func- 
tional associated with the particular shape is related to 
the free energy determined operationally. 

[13] If there is no conservation constraint, the combined prob- 
ability of Xk and X n ~k, P(Xk,X n -k), approaches a 
product of Gaussian distributions as k and n — k be- 
come large. The asymptotic distribution of Xk with the 
conservation constraint is simply given by P(Xk, —Xk)- 

[14] The steady state entropy operationally defined by Oono 
and Paniconi |7j is related to the Shannon entropy. See 
T. Hatano and S. Sasa, Phys. Rev. Lett. 86, 3463, (2001). 



